module mod_elglike_element
use mod_matrix
!use mod_shearpower
  implicit none
 
  !
  ! Elg-likelihood
  !
  ! Contains methods for conversion between various structures, convolusion, automatic loading and saving of maps & power spectra,
  ! timer functions, parameter file usage, output and misc
  !
  ! (c) Nicolaas E. Groeneboom
  !     leuat@irio.co.uk
  !     http://www.irio.co.uk
  !

  type, public :: t_elglike_element
     type(t_matrix) :: array
     real*8 :: determinant
     real*8 :: likelihood
  end type t_elglike_element
 
  integer :: zbins = 2


double precision :: global_KN = 0.22**2 *3600*180/(3.14159**2)/0.9*1.412

 
contains  

subroutine elglike_element_initialize(elglike_element)
  type(t_elglike_element), intent(inout) :: elglike_element
  call matrix_initialize(elglike_element%array, zbins)
  elglike_element%likelihood = 0
  elglike_element%determinant = 0
end subroutine elglike_element_initialize

subroutine elglike_element_deallocate(elglike_element)
  type(t_elglike_element), intent(inout) :: elglike_element
  call matrix_deallocate(elglike_element%array)
  elglike_element%likelihood = 0
end subroutine elglike_element_deallocate



function elglike_element_zfunc(z)
  double precision :: elglike_element_zfunc
  double precision, intent(in) :: z

  elglike_element_zfunc = exp(-z);

end function elglike_element_zfunc

!!$function elglike_element_calculate_N()
!!$  double precision :: elglike_element_calculate_N, val,a,b
!!$  integer :: i,j
!!$
!!$  val = 0
!!$  do i=1, nbins-1
!!$     a = zbins(i)
!!$     b = zbins(i+1)
!!$
!!$     val = val + (b-a)*0.5*(elglike_element_zfunc(a) + elglike_element_zfunc(b))
!!$  end do
!!$  elglike_element_calculate_N = global_KN / val
!!$  
!!$
!!$end function elglike_element_calculate_N 

function elglike_element_likelihood(elg, invmat, data, l)
  real*8 :: elglike_element_likelihood
  real*8, intent(in) :: l
  real*8 :: trace
  type(t_matrix), intent(inout) :: invmat
  type(t_elglike_element), intent(inout) :: elg, data
  type(t_matrix)  :: tempmat
  

  call matrix_initialize(tempmat,zbins)

  
  tempmat%val = matmul(invmat%val, data%array%val)
  trace = matrix_trace(tempmat)

  trace = trace + log(elg%determinant/data%determinant)
  trace = trace - zbins

  trace = trace * (2.0*l+1)

  call matrix_deallocate(tempmat)
  
  elglike_element_likelihood = trace

end function elglike_element_likelihood


subroutine elglike_element_setup_matrix(elg, param)
  real*8 :: l, tmp
  integer :: i,j
  double precision, intent(in) :: param(10)
  type(t_elglike_element), intent(inout) :: elg
  type(t_matrix) :: invmat, testmat
  real*8 :: d


!  d = elglike_element_calculate_N()
!  write(*,*) "N: ", d
  d = 5.85E-10 
  
  l = param(6)
  !write(*,*) "Setting up for ", param
  do i=1, zbins
     do j=1,i
        !write(*,*) i,j,l, elg%array%val(i,j)
        call shearpower(param(1:5), i, j,l, elg%array%val(i,j))

        if (i==j) elg%array%val(i,j) = elg%array%val(i,j)  + d
        elg%array%val(j,i) = elg%array%val(i,j)
     end do
  end do

  elg%array%val = elg%array%val * 1E8
end subroutine elglike_element_setup_matrix


subroutine elglike_element_calculate_determinant(elg, l)
  real*8, intent(in) :: l
  type(t_elglike_element), intent(inout) :: elg
  type(t_matrix) :: invmat
  
  call matrix_initialize(invmat, zbins)
  call matrix_invert(elg%array, invmat)
  elg%determinant = invmat%determinant

end subroutine elglike_element_calculate_determinant

subroutine elglike_element_calculate_likelihood(elg, l, data)
  real*8, intent(in) :: l
  integer :: ll
  type(t_elglike_element), intent(inout) :: elg
  type(t_elglike_element), intent(inout) :: data
  type(t_matrix) :: invmat, testmat
  
  call matrix_initialize(invmat, zbins)

  call matrix_invert(elg%array, invmat)

  elg%determinant = invmat%determinant

  ! need to be data!

  ! change to DATA
  ll = l
  elg%likelihood = elglike_element_likelihood(elg, invmat, data, l)



 ! call matrix_initialize(testmat, N)
 ! testmat%val = matmul(orgmat%val, invmat%val)
 ! write(*,*) "multiplied:"
 ! call matrix_print(testmat)
 ! write(*,*) "orgmat:"
 ! call matrix_print(elg%array)

  
  call matrix_deallocate(invmat)

end subroutine elglike_element_calculate_likelihood




end module mod_elglike_element
